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Abstract 

Empirical studies over the past two decades have provided support for the hypothesis that 
schizophrenia is characterized by altered connectivity patterns in functional brain net- 
works. These alterations have been proposed as genetically-mediated diagnostic biomark- 
ers and are thought to underlie altered cognitive functions such as working memory. 
However, the nature of this dy connectivity remains far from understood. In this study, 
we perform an extensive analysis of functional connectivity patterns extracted from MEG 
data in 14 subjects with schizophrenia and 14 healthy controls during a 2-back working 
memory task. We investigate uni-, bi- and multivariate properties of sensor time series by 
computing wavelet entropy of and correlation between time series, and by constructing 
binary networks of functional connectivity both within and between classical frequency 
bands (7, /?, a, and 9). Networks are based on the mutual information between wavelet 
time series, and estimated for each trial window separately, enabling us to consider both 
network topology and network dynamics. We observed significant decreases in entropy in 
prefrontal and lateral sensor time series and significant increases in connectivity strength 
in the schizophrenia group in comparison to the healthy controls. We identified an in- 
verse relationship between entropy and strength across both subjects and sensors that 
varied over frequency bands and was more pronounced in controls than in patients. The 
topological organization of connectivity was altered in schizophrenia specifically in high 
frequency 7 and j3 band networks as well as in the 7-/3 cross-frequency networks. Network 
topology varied over trials to a greater extent in patients than in controls, suggesting 
disease-associated alterations in dynamic network properties of brain function. Our re- 
sults identify signatures of aberrant neurophysiological behavior in schizophrenia across 
uni-, bi- and multivariate scales and suggest novel candidate intermediate phenotypes in 
cross-frequency network architecture and in network dynamics. 
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Introduction 

Synchronous oscillatory brain activity is thought to form a neurophysiological net- 
work over which the brain processes information |41| during cognitive processes as varied 
as learning and memory. Coherent neuronal communication is constrained by the un- 
derlying neurochemistry of the dopamine, GABA, glutamate and acetylcholine systems 
[2"2l I51] . Diseases like schizophrenia thought to involve such neurotransmitter systems 
|22l[32l[94] are therefore also often accompanied by altered oscillatory network structure. 

Indeed, a hallmark of schizophrenia is the complex pattern of abnormal increases and 
decreases in resting state and task-based connectivity evident between large-scale brain 
regions EU ESI |H| . To quantify this dysconnectivity [35] at the level of the whole brain, 
network theoretical tools originally developed in the social sciences have been applied to 
neuroimaging data. Results demonstrate complex dysconnectivity profiles of the brain's 
anatomy [HH1 Ell EE] and function [TQ1 H3j Ml Ml EH Ml EH Ell E] in people with 
schizophrenia [55] , 

Altered brain function in schizophrenia is unsurprisingly coupled with altered cog- 
nitive behavior. Due to the heterogeneity of the disease, some symptoms including 
both positive (hallucinations, delusions and thought disorder) and negative (e.g. im- 
paired emotional response and social interaction, avolition, anhedonia) symptoms, vary 
strongly among patients. However, cognitive impairments in executive functions (such 
as the manipulation of transiently stored information [55 ) and working memory (e.g., 
the ability to transiently maintain and manipulate a limited amount of information in 
order to guide thought or behavior) are typically altered to varying degrees in virtually 
all patients. These cognitive impairments appear early and change little over time, and 
are therefore thought to lie at the core of the disease, contributing to the development 
of clinical symptoms [331 144) . Both altered memory and executive planning have been 
linked to the neurochemistry of dopamine, GABA, glutamate and acetylcholine systems 
and subsequently to aberrant (NMDAR)-mediated synaptic plasticity [351 12H EI] • 

Correlations between brain network organization and cognitive or behavioral variables 
in schizophrenia are expected |28j but quantitative evidence for such a relationship is 
still rare. Two recent studies have demonstrated that the organization of task-based 
functional brain networks in this population can be linked both to the efficacy of auditory 
rehabilitation efforts EE] and to working memory performance [T31 051 EH E3] ■ These 
findings suggest that the characterization of the dynamic network of synchronous brain 
activity might be a powerful tool to study cognitive impairments in schizophrenia and 
might provide insight into the neurobiological mechanisms of the disease. However, 
many of these studies have several potential limitations, such as using only a small set 
of network diagnostics and constructing a single functional network from a large array 
of data. 

In this study, we provide an expanded examination of the functional brain network ar- 
chitecture in people with schizophrenia estimated from magnetoencephalography (MEG) 
data acquired during the performance of a working memory N-back task. We first exam- 
ine the link between single time series variability and pairwise time series similarity, as 
a means of characterizing the relationship between functional activity and connectivity. 
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To quantitatively characterize functional brain network structure, we examine 12 binary 
network diagnostics across a wide range of densities and employ a recently developed 
statistical technique (functional data analysis) to identify significant group differences. 
To detect time-varying network structure and to increase statistical robustness, we con- 
struct an ensemble of networks for each individual using data extracted from 66 separate 
trials. 

Finally, in a novel extension of the analysis of functional brain network structure in the 
classical frequency bands of brain activity, we examine functional networks constructed 
from the interactions between frequency bands. Inter-frequency networks are particularly 
interesting in the context of our examination of the N-back task because cross-frequency 
interactions are thought to facilitate memory function |79j by enabling the cross-function 
integration of information over spectrally distinct processing streams [55]. Furthermore, 
evidence suggests that cross-frequency interactions are altered in schizophrenia, and the 
strength of these alterations has been linked to genetic risk factors for the disease [5J. 

Using these tools, we uncover an extensive pattern of altered network structure and 
network dynamics in people with schizophrenia that could potentially be used to identify 
intermediate phenotypes and develop diagnostic biomarkers for the disease. 

Materials and Methods 

Data Acquisition 

14 healthy volunteers and 14 people diagnosed with schizophrenia spectrum disorders 
(according to the Diagnostic and Statistical Manual of Mental Disorders IV criteria) took 
part in the Clinical Brain Disorders Branch/National Institute of Mental Health Genetic 
Study of Schizophrenia (National Institutes of Health Study Grant NCT 00001486, Daniel 
R. Weinberger, principal investigator). None of the healthy volunteers had structural 
magnetic resonance imaging abnormalities or history of psychiatric illness, depression, or 
loss of consciousness. All patients were receiving antipsychotic medication at the time of 
the study; none of the healthy volunteers were taking psychoactive medication. Subjects 
were matched for sex (9 males, 5 females in each group) and did not differ significantly in 
age. The age of the healthy control group was 30 years ±7.3 (SD), and that of the patient 
group was 33 years ±8.6 (SD). The protocol was approved by the National Institute of 
Mental Health Institutional Review Board and informed consent was given in writing by 
all participants. 

MEG data were acquired at the National Institute of Mental Health in Bethesda, MD 
using a 275-channel CTF system (VSM MedTech) with a sampling rate of 600 Hz. The 
experimental paradigm was a visual 2-back working-memory task using numerical stimuli 
ranging from 1 to 4 as previously described |24j . Earlier studies using the same task 
demonstrated abnormal behavioral and neuroimaging associations with schizophrenia 
[H [23l [67]. Each subject performed six blocks of 11 trials (66 trials in total per subject). 
In each trial, a number was presented visually for 500 msec and, beginning with the third 
trial, the subject was asked to respond within 1,300 ms by pressing 1 of 4 buttons to 
indicate the identity of the number seen 2 trials previously. 

Data Preprocessing 

Using MATLAB [58 and FieldTrip software [51] , raw data were mean corrected and 
filtered to attenuate background low-frequency noise and line noise at 60 Hz by using 
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a 0.3-Hz-width filter. The axial gradiometers of the CTF machine have source profiles 
that include information from a wide spatial range, strictly limiting the interpretation 
of the anatomical location of results. To retain the inherent correlation structure of a 
network of interacting brain regions while gaining localization specificity, we transformed 
the data into planar space. Using Fieldtrip, the planar transform was applied to the data 
by using the function ft.megplanar. For trial-by-trial analysis, the time series were cut 
into epochs of 1.8s, beginning with the presentation of one stimulus and ending with the 
beginning of the next stimulus (or with the end of the block of trials in the case of the 
final trial). We therefore investigated a total of 66 time windows per subject. 

Wavelet Decomposition 

We focused our investigation on frequency-band-specific oscillations in the MEG sig- 
nal by employing a wavelet analysis. Time series were resampled to 120 Hz to constrain 
the frequency bands of the wavelet transform to roughly conform to the classical fre- 
quency bands [13] . Each frequency band of interest was isolated by applying the maximal 
overlap discrete wavelet transform to each time series [66] . In line with previous work, 
we used the Daubechies 4 wavelet [13]. Wavelet scale 1 (30-60 Hz) roughly corresponds 
to the 7 band, scale 2 (15-30 Hz) to j3, scale 3 (8-15 Hz) to a, and scale 4 (4-8 Hz) to 9. 
Because of the relatively short trial length (1.8 s), we did not examine frequencies below 
4 Hz. 

Functional Connectivity and Complexity Estimation 

Before constructing functional brain networks, we characterized each sensor's wavelet 
time series and the pairwise relationship between sensor time series. To quantify the 
complexity of the time series, we computed the Shannon wavelet entropy [12 [731 E] > 
which we describe in detail in the supplementary online materials (SoM). To investigate 
functional connectivity between sensors, we calculated the mutual information (MI), 
which is particularly appropriate for estimating interactions that encompass a narrow 
frequency band (e.g., a wavelet scale) [25] . 

We represent pairwise sensor MI both within and between frequency bands as func- 
tional connectivity matrices. To construct intra-frequency functional connectivity ma- 
trices, we calculated the MI between the time series in frequency band a of all possible 
pairs of sensors to create the weight matrix W a . To construct inter- frequency functional 
connectivity matrices, the process was similar: for each pair of frequency bands a and b, 
we calculated the MI between the time series of sensor i in frequency a and sensor j in 
frequency b for all possible pairs of sensors i and j. This resulted in the weight matrix 
W ab . The pairwise mutual information values in both the intra- and inter-frequency 
matrices were normalized according to Strehl and Ghosh |84j , ensuring that the values of 
the elements of W a and W ab were in the range [0, 1]. Through this process, we created 
4 intra-frequency functional connectivity matrices (7, /3, a, and 9) and 6 inter-frequency 
functional connectivity matrices (7 — ft, 7 — a, 7 — 9, /? — a, P — 0, a — 9) per subject for 
each of the 66 trials (660 total networks per subject). 

The strength of sensor i was defined as the average MI between that sensor and all 
other sensors: Si = ^ Y^j The strength of a participant was defined as the average 

strength of all sensors: S = ^ Si. 



4 



Network Properties 

To quantitatively characterize the topology of the intra- and inter-frequency func- 
tional connectivity matrices, we construct sets of binary networks summarizing the in- 
teractions (edges) between sensors (nodes). A binary network can be represented math- 
ematically by an adjacency matrix A, whose entries Aij are either or 1, indicating the 
absence or presence of an edge respectively. Binary networks can be obtained by thresh- 
olding a given functional connectivity matrix in several ways |81l 114] . Here we employ 
the technique of cumulative thresholding, a procedure in which a threshold is applied 
to the weighted functional connectivity matrix W to retain a given percent of strongest 
connections. Weights (e.g., MI values) that pass this threshold were set to a value of 1 in 
the adjacency matrix A, while those that did not pass this threshold were set to a value 
of in the adjacency matrix A. 

The percent of nonzero elements of A is also known as the network density, or network 
cost K.. We can employ a range of thresholds to probe the topological organization of a 
network over a range of network densities. Using a high threshold, we obtain a sparse 
matrix representing the topology of the strongest functional associations between sensors, 
while using a lower threshold, we obtain a denser matrix representing the topology of 
functional associations over a wider range of strengths. In this study, we employed a 
set of thresholds enabling us to examine corresponding sets of networks with densities 
ranging from (no connections present) to 0.5 (half of the possible number of connections 
present) in steps of 0.01. We refer to this network density range < JC < 0.5 as the cost 
regime of interest. Our choice to examine network organization over a large range of 
network density values is supported by previous studies EE] that have demonstrated 
that the choice of threshold can have a large influence on the topological properties of 
the binary graph. 

We quantified the organization of the binary networks using the path-length and 
clustering coefficient [95] , global and local efficiency [531 H US] , betweenness centrality 
|gD] . modularity [HIE!], hierarchy [THOU], synchronizability 0Q1], assortativity [551 ITU] , 
and robustness [2j [57] . In addition to topological network properties, we also studied the 
physical measure of mean connection distance |10j . as well as a combined topophysical 
property Rent's exponent [TTJ which estimates the efficiency of the topological embedding 
of the network into physical space. Mathematical descriptions of these diagnostics are 
given in the SoM. 

Statistical Analyses 

We examined the reliability of network properties over subjects and trials using the 
coefficient of variation (CV), which is defined as the standard deviation a of a given 
sample normalized by its mean value \x: 

CV=-. (1) 
A 4 

To test for inter-subject temporal variability in functional brain network organization, 
we calculated the CV for each network diagnostic, for each subject in both groups, and 
for all 10 intra- and inter-frequency networks. Note that we excluded the hierarchy and 
assortativity from this analysis because their values are close to zero, making the CV less 
meaningful and more prone to estimation errors from the division of numbers < 1. To 
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test for group differences in temporal variability of network structure, we used a repeated 
measures ANOVA with CV values averaged over costs, group as a categorical factor, and 
frequency band and graph diagnostic as repeated measures. 

To identify group differences in network diagnostics, we used Functional Data Analysis 
(FDA). FDA enables statistical inference from sets of functions [70] by extending the 
principles of statistical inference from data points to data curves. Here, in keeping with 
|14j . the values of the graph properties were treated as a function of network density, 
and the two groups, (controls and people with schizophrenia) were compared with a 
non-parametric permutation test using twenty thousand permutations of group labels. 

In addition to testing for group differences in network structure, it is often of interest 
to determine whether the structure of an empirical network is different from what one 
would expect in a given null model. While the development of potentially useful network 
null models is ongoing [1031 115j , here we employ benchmark Erdos-Renyi (ER) random 
graphs to test whether the network topology identified in intra- and inter-frequency 
functional brain networks was non-random. We created an ensemble of 66 ER graphs 
for each network density and each of the 10 frequency bands. Note that the number of 
networks in the ensemble was set to be identical to the number of networks in a single 
subject. 

All computational and statistical operations were implemented in MATLAB® (2007a, 
The MathWorks Inc., Natick, MA). Network diagnostics were estimated using a com- 
bination of in- house software, the Brain Connectivity Toolbox [76], and the MATLAB 
Boost Graph Library [58]. The repeated measures ANOVA was performed using freely 
available code |106j . 

Results 

Working Memory Performance 

We examined the accuracy of performance in the 2-back working memory task for 
both the schizophrenia and control groups. People with schizophrenia had a significantly 
lower accuracy (44.4±25.5% (STD), the median was 56.1%) than controls (84.2 ± 16.7%): 
two sample t-test t = 5.15, p = 2.23 x 10~ 5 . These results confirm that working memory 
function is impaired in our cohort of people with schizophrenia, supporting an additional 
investigation into the patterns of brain function during task performance. 

Time Series Variability and Co- Variability 

Time Series Variability. We examined the complexity of the MEG signals in the two 
groups by measuring the Shannon wavelet entropy of the sensor time series in four clas- 
sical frequency bands(7-, /?-, a- and 8; see Methods). On average, the entropy was lower 
in people with schizophrenia in 3 of the 4 bands (7, (3, and a); see Figure [TJ\_. Further- 
more, the distribution of entropy across sensors varied over frequency bands and between 
groups; see Figure [2] The most significant group differences were identified in lateral and 
prefrontal areas, where entropy is lower in people with schizophrenia than in healthy 
controls in the high frequency 7- and /?- bands. 
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Figure 1: Group Differences in Activity and Connectivity. Wavelet entropy (A) and intra- 
frequency strength (B), averaged over sensors, in the 7-, f3-, a- and 0-bands for healthy controls and 
people with schizophrenia. Asterisks indicate significant group differences as measured by two-sample 
t-tests (p < 0.05). 




Figure 2: Anatomical Distribution of Entropy. Topography of wavelet entropy estimated for 
individual sensors in the 7-, /3-, a- and 6-bands (A-D respectively), averaged over the people with 
schizophrenia (Top row) and healthy controls (Second row). Difference between groups (probands minus 
controls) (Third row) were tested for significance (p < 0.05) by permutation testing of group labels 
(Bottom row). 
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Figure 3: Anatomical Distribution of Strength. Topography of average mutual information between 
individual sensors and all other sensors in the 7-, /?-, a- and 0-bands (A-D respectively), averaged over 
the people with schizophrenia (Top row) and healthy controls (Second row). Difference between groups 
(probands minus controls) (Third row) were tested for significance (p < 0.05) by permutation testing 
(Bottom row). 



Time Series Co-Variability. The co- variability of sensor time series is thought to be 
a measurement of synchronous oscillatory neuronal activity and therefore potentially a 
proxy for large-scale communication between brain regions. We estimated the mutual 
information between sensor time series in the same four frequency bands (7-, /3-, a- and 
9; see Methods). The strength of a sensor was defined as the average mutual information 
between that sensor's time series and the time series of all other sensors. Strength 
displayed frequency-dependent variation, being smallest in the highest frequency band 
(7) and largest in the lowest frequency band (9) for both groups; see Figure [T] On 
average, the strength was higher in people with schizophrenia in 3 of the 4 bands (7, a, 
and 9). Furthermore, the distribution of strength across sensors varied over frequency 
bands and between groups; see Figure [3j We found that strength was larger in people 
with schizophrenia, particularly in prefrontal sensors in the 7-band and in prefrontal and 
other lateral sensors in the f3- band. 

Relationship Between Entropy and Strength. Recent work has suggested a potential and 
intuitive link between time series variability (e.g., entropy) and time series co- variability 
(e.g., correlation or mutual information) in fMRI data over subjects 104 and over brain 
regions |14| . Here we examine this potential relationship in a different imaging modality 
(MEG) both across frequency bands and between groups. In Figure |9j we show the 
strength and entropy of each individual, averaged over sensors and trials. Importantly, 
we find that correlations between the two measurements at the inter-subject level are 
frequency-dependent. Both groups display a significant inverse relationship between 
entropy and strength in the /3-band, suggesting that individuals with high temporal 
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Figure 4: Entropy and Strength: Inter-Subject Level. Scatterplots of the strength of connectivity 
and the complexity of the time series, as measured by wavelet entropy, for 7-, /3-, a- and 0-bands. 
Single data points represent values for each individual averaged over trials and sensors. Red markers 
denote subjects with schizophrenia spectrum diagnosis; blue markers denote healthy subjects. Lines 
indicate best linear fits for the two groups separately (red, SZ; and blue, HC) and we provide r 2 values 
as indicators of goodness of fit. Similar scatterplots that code each experimental block separately are 
provided in Figure SI in the SoM. 



variability in brain function have lower temporal co-variability. This relationship also 
appears to hold for the controls in the other frequency bands. In Figure [5] we show the 
strength and entropy for each sensor (rather than each individual), averaged over trials 
and over individuals within a group. Again, we identify a strong inverse relationship 
between entropy and strength, which in this case is strongest in the 0-band. Of note, 
the inhomogeneous distribution of values characterized by a sparse low entropy tail is 
consistent with data reported in [14], and potentially indicates variation in the roles of 
brain areas in information processing. 

Network Structure 

Given the group differences in connectivity strength identified in the previous sec- 
tion, we next asked whether the patterns of connectivity between sensors are altered 
in schizophrenia. We characterize these patterns using binary network diagnostics. See 
Methods for details on binary network construction and the SoM for mathematical defi- 
nitions of network diagnostics. 

Cost-Efficiency. We constructed binary graph diagnostics as a function of network den- 
sity using a cumulative thresholding technique (see Methods). A simple network diag- 
nostic that collapses such a curve into a single value is the cost-efficiency, which is defined 
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Figure 5: Entropy and Strength: Sensor Level. Scatterplots of the strength and wavelet entropy 
in 7-, j3-, a- and 0-bands. Single data points represent values for each sensor averaged over trials and 
individuals. Red markers denote data from the SZ group; blue markers denote data from the HC group. 
Lines indicate linear fits for the two groups separately (red, SZ; and blue, HC) and we provide r 2 values 
as indicators of goodness of fit. 
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Figure 6: Cost-Efficiency of Functional Networks in Health and Disease. (A) Cost-efficiency 
in the 7-, /?-, a-, 8-, and cross- frequency band networks for healthy controls (blue) and people with 
schizophrenia (red). Error bars indicate the standard error over subjects. The gray shaded line indicates 
the expected values of cost-efficiency for an ensemble of random (ER) graphs (see Methods). 



as the maximum of the efficiency-minus-cost curve [TJ [T31 EH] (see SoM). In Figure [6]A., 
we show the cost-efficiency of both intra- and inter-frequency networks for both groups. 
As expected from previous work [T3] , we found that cost-efficiency values decreased with 
increasing frequency. Importantly, intra-frequency networks demonstrated consistently 
lower cost-efficiency than inter- frequency networks. Group differences were only evident 
in intra-frequency networks (e.g., the /3, a and 9 networks), with cost-efficiency being 
higher in the control group than in the patient group. Inter-frequency networks, while 
not demonstrating a group difference, did show a significantly higher cost-efficiency than 
expected in an ensemble of Erdos-Renyi random graphs (see Methods), suggesting the 
presence of non-random structure in cross-frequency interactions. 

Binary Network Organization. While cost-efficiency has the advantage of collapsing a 
binary-diagnostic versus cost curve into a single value, it is also of interest to examine 
the shape of these curves for the other diagnostics. We examined 12 binary graph diag- 
nostics as a function of cost for both groups and all 10 frequency bands [8]; see Methods 
and SoM. We separated diagnostics into those that showed higher values in people with 
schizophrenia (e.g., see Figure [7|\) and those that showed higher values in healthy con- 
trols (e.g., see Figure [7j3) . We found that the majority of significant group differences 
were located in the 7 and f3 intra-frequency networks and in the 7-/? cross-frequency 
networks; see Figure [7p. Note that the a-band also showed group differences for 5 out 
of the 12 diagnostics. People with schizophrenia displayed higher global and local effi- 
ciency, betweenness centrality, clustering coefficient, modularity, and assortativity and 
displayed lower hierarchy, synchronizability, mean connection distance and robustness 
to both targeted and random attack. The identification of group differences in both 
topological and physical network properties suggests that disease-associated changes in 
brain network organization might be linked to multiple developmental mechanisms. In- 
deed, recent theoretical work has suggested that altered constraints on both information 
efficiency (a topological property) and metabolic cost (a physical property) can lead to 
schizophrenia- like changes in network architectures |89) . 
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Figure 7: Binary Network Organization in Health and Disease (A-B) Example binary net- 
work diagnostic curves as a function of threshold: the modularity for controls (black) and people with 
schizophrenia (red) (A), and the robustness to targeted attack for controls (black) and people with 
schizophrenia (blue) for 7-band networks. Individual curves indicate average values for each individual 
over the 66 trial-specific networks. (C) Significant group differences in graph diagnostic versus cost 
curves for 12 graph diagnostics (y-axis) in both intra- and inter-frequency bands (x-axis). Warm colors 
indicate that the diagnostics values were higher in people with schizophrenia than in healthy controls 
(red, p < 0.01; orange, p < 0.05 uncorrected); cool colors indicate that the values were higher in healthy 
controls than in people with schizophrenia (dark blue, p < 0.01; light blue, p < 0.05 uncorrected). See 
SoM for the full graph diagnostic versus cost curves for all subjects, frequency bands, and diagnostics. 
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Figure 8: Temporal Variability of Network Diagnostics for binary network diagnostics in the 7 
and /9 intra-frequency networks and in the 7-/3 inter-frequency networks. CV values indicate temporal 
variability over trials for healthy (blue) and schizophrenic subjects (red), averaged over the entire range 
of cost values. Error bars indicate the mean squared error over subjects and costs. Values and error 
bars for synchronizability are scaled down by a factor of 10 for visualization purposes. Results for other 
frequency bands can be found in the SoM. 



Variability of Binary Network Structure. Thus far we have reported results for each 
individual derived from networks constructed from 66 trial blocks. Here we ask whether 
network organization is variable over trial blocks in the two groups. To quantify temporal 
variability in network structure, we computed the coefficient of variation (CV) for network 
diagnostics over trial blocks for each individual in all 10 frequency bands. 

In Figure [12] we show the CV for the 7 and (3 intra-frequency and the 7-/3 inter- 
frequency networks, averaged over costs (for similar results in other frequency bands, see 
the SoM). We observe that the CV varies over binary graph diagnostics in a similar way 
across frequency bands. We also note that in almost all cases, the temporal variability 
appears to be larger for people with schizophrenia. To quantitatively test this observa- 
tion, we used a Repeated Measures ANOVA (see Methods for details and Table 1 for 
results). The main effect of group (F = 6.24, p — 0.02) confirmed that the networks 
derived from people with schizophrenia varied more over time than did those derived 
from healthy controls, suggesting a fundamental alteration in the dynamics of functional 
brain networks in schizophrenia. 

Discussion 

We have examined temporal characteristics of MEG data acquired from people with 
schizophrenia and controls during a working memory task. Our approach spanned sev- 
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Table 1: Results for a Repeated Measures ANOVA of Temporal Variability in Network 
Topology, measured using the CV. 
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eral distinct levels, including that of the individual sensor time series (univariate), the 
co- variability between time series (bivariate), and the patterns of co- variability across 
sensors using binary network diagnostics (multivariate) . We identified disease-associated 
changes in brain function at each level. People with schizophrenia displayed lower time 
series entropy, higher strength of co- variability between time series, and an extensive pat- 
tern of altered topological organization in binary sensor networks. Importantly, network 
properties of cross-frequency associations between time series in the 7- and /3-bands dif- 
fered between groups, uncovering a novel feature of dysconnectivity. Finally, the temporal 
variability of brain network architecture in people with schizophrenia was significantly 
higher than that in healthy controls, a phenomenon suggestive of decreased dynamic 
constraints on brain function. 

Identification of Intermediate Phenotypes 

The mechanisms of psychotic symptoms in schizophrenia, their genetic underpinnings 
and their brain signatures, are far from understood. The interpretations of neuroimag- 
ing phenomena and related conclusions regarding neurophysiological mechanisms of the 
disease, are greatly hampered by the confounding effects of medication and disease het- 
erogeneity associated with epiphenomena related to smoking and chronicity, as well as the 
heterogeneity of the mental state of ill subjects during these MEG procedures. Thus, it 
is impossible to conclude that our findings represent primary disease phenomena rather 
than epiphenomena related to secondary factors that are associated with the state of 
illness. These difficulties can be overcome to some degree by leveraging the genetic sim- 
ilarity between people with schizophrenia and their unaffected siblings and testing for 
similar imaging phenomena in both groups. 

An important goal of this study was to perform an extensive assessment of alter- 
ations in brain functional architecture in schizophrenia to identify candidate diagnostics 
for intermediate phenotypes of the disease |36| . Intermediate phenotypes - observable 
characteristics that show tiered values in patients, their siblings, and controls - represent 
a uniquely powerful tool for the quantification of genetic risk mechanisms for psychiatric 
disease [60 . In particular, phenotypes derived from human neuroimaging have intriguing 
potential to provide concrete insight into the neurochemical pathologies |47j that might 
directly underlie the behavioral symptoms of psychosis |86j . 
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Our multiscale exploration of illness-associated physiologic phenomena and their po- 
tential to represent intermediate phenotypes has been guided by evidence for genetic 
regulation of both brain activity variability [97] and functional connectivity architecture 
|37j . Temporal variability of brain activity specifically during memory tasks has previ- 
ously been identified as an intermediate phenotype of schizophrenia |97) . validating the 
use of the iV-back task in our study. Further examining the patterns of functional con- 
nectivity supporting memory function is critical given recent evidence suggesting that 
functional connectivity can be more sensitive to genetics [M] and to schizophrenia [T^l[T5] 
than the simple activity of brain regions alone. 

Our results suggest that network organization in the higher frequency bands (7 and 
f$) is potentially a better candidate for an intermediate phenotype than that in lower fre- 
quency bands during the ./V-back task. Importantly, in this study, we have also introduced 
two other potential candidates for intermediate neuroimaging phenotypes based on sig- 
natures of dynamic sensor-sensor communication: 1) network properties of whole-brain 
inter-frequency connectivity and 2) dynamic changes in network configuration during task 
performance. We find that the 7-/3 networks display significant group differences in orga- 
nization, similar to the patterns seen in the high frequency bands alone. Furthermore, the 
variability in this network structure over the duration of task performance is significantly 
higher in people with schizophrenia than in controls, suggesting that dynamic reconfig- 
uration properties of the healthy brain might be affected in the disease. Taken together, 
our results lay the groundwork for a directed examination of the genetic underpinnings 
of synchronous oscillatory brain activity in probands and their unaffected siblings. Im- 
portantly, as deficits in cognitive function such as working memory are thought to be 
fundamental and possibly causal in schizophrenia 33, 13] , such studies might eventually 
lead to improved detection and treatment of the disease. 

Altered Time Series Variability 

Signal variability or noise is a characteristic feature of the cortical system and 
is thought to facilitate the exploration of functional network configurations necessary 
for healthy cognitive function |31j . Increased variability of stimulus- induced prefrontal 
electromagnetic activity has been observed in both people with schizophrenia and their 
healthy siblings |101l Infill 00) ■ suggesting that cortical noise might be a genetic biomarker 
for the disease. Indeed, distributed patterns of both increased and decreased signal 
variability differentiate people with predominantly positive or negative symptom profiles 
|69| 140] . Furthermore, signal variability is behaviorally relevant, having been linked to 
task accuracy in both probands and healthy controls H7J 1521 I5T)] . 

Our results support the notion that signal variability is affected in schizophrenia. By 
employing the simple measure of time series entropy, we showed that signal variability 
was decreased in probands specifically in prefrontal and lateral sensors in high frequency 
bands (7 and j3). However, it is important to note that the Shannon entropy character- 
izes properties of the distribution of values in the time series rather than the temporal 
evolution of the signal. Alternative measures of temporal signal variability, such as the 
multi-scale entropy |85| , the Renyi number [45] , the Lyapunov exponent [102] , the fractal 
dimension (77JH25], and the Hurst exponent [2D], might provide greater or lesser sensitivity 
to disease state [78] or to disease-associated alterations in lower frequency bands. 
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Activity and Connectivity 

Intuitively, one might imagine that the complexity of a region's time series would have 
some bearing on the strength of connectivity between that region and other regions. For 
example, a region with a more noisy signal might be less likely to show strong associations 
with other regions than a region with an ordered, strong oscillatory signal. 

Two recent fMRI studies have supported the view of a strong positive correlation 
between time series variability and co- variability [LH 1104] . However, our results suggest 
a more complex relationship in MEG data that is modulated by both frequency and 
disease. At the subject level, an inverse correlation between entropy and strength was 
evident in the (3 band for both groups, and also in the 7, a and 9 bands for the control 
group. Similarly at the sensor level, the entropy and strength were again inversely corre- 
lated, but in this case the effect was largest in the 9 band. Importantly, the sensor-level 
analysis uncovered a strong heterogeneity in the distribution of entropy across sensors, 
complicating the statistical estimate of the entropy-strength relationship. 

The consistent inverse relationship between entropy and strength identified in this 
study is at odds with the positive relationship identified in previous fMRI studies [T41I104) . 
In our case, high entropy sensors tend to have weak connectivity to other sensors, while 
in the two previous fMRI studies, regions or individuals characterized by high entropy 
tended to also show strong connectivity. However, the divergence between results from 
MEG and fMRI studies might not be that surprising, as the two imaging modalities 
measure inherently different properties of neurophysiological activity (blood flow and 
magnetic flux) at different temporal resolutions. The two imaging modalities also appear 
to capture signatures of brain communication differently; for example, increases in the 
BOLD signal have been linked to both synchronization and desynchronization in MEG 
oscillations |99j . Further work in both imaging modalities will be necessary to gain a 
greater intuition for the role of activity-connectivity relationships in cognitive function. 

Network Structure of Connectivity 

Both binary and weighted network analyses each have their specific advantages and 
disadvantages and, when applied to the same data, can yield diverging results [151 126j . 
Weighted networks suffer from issues of normalization, which are particularly important 
in the context of comparing groups with different average strength (as is the case here in 
our study), and further require the careful identification of appropriate null models [761 
103, 15j. Binary networks on the other hand rely heavily on the choice of thresholds, and 
neglect information contained in the weights of connections. However, binary networks 
are mathematically simpler |76j and their disadvantages can be at least partially overcome 
by examining diagnostic properties over a range of network densities. In this paper, 
we focus on binary networks and examine their properties over a wide and dense cost 
range. To identify group differences, we employ functional data analysis |14| and find 
that people with schizophrenia display higher local efficiency, betweenness centrality, 
clustering coefficient, modularity, and assortativity and display lower global efficiency, 
hierarchy, synchronizability, mean connection distance and robustness to both targeted 
and random attack than healthy controls in high frequency band networks. 

In comparing our findings to the literature, we note that results of connectivity anal- 
yses seem to differ between imaging modalities. Network organization in schizophrenia 
has been examined using a variety of structural (sMRI, DTI [TUl 11041 19"3"] ) and functional 
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(fMRI, EEG, and MEG) neuroimaging techniques [TH HH EU E21 [571 HS1 El II]- While 
schizophrenia-related alterations in network structure are reported across all imaging 
modalities and experimental paradigms, the patterns of these differences do not always 
converge. One relatively consistent finding is that global network efficiency is lower 
in schizophrenia subjects for networks derived from DTI 104, 93 , resting state fMRI 
data [56l [92] , and our current MEG study. However, the clustering coefficient has been 
reported to be lower in schizophrenia for several resting state fMRI and EEG studies 
[751 [571 [55] but is higher in our investigation of task-based MEG in one inter- and several 
intra-frequency networks. While this divergence of results could be due to the different 
imaging modalities, another possible interpretation might be that the direction of net- 
work organizational changes in schizophrenia could depend on brain state (e.g., rest or 
task) . Such a view is supported by recent work demonstrating that network organization 
changes as a function of cognitive load [J51 [5T] and learning [T5] . 

The complex landscape of network changes in schizophrenia demonstrates the need 
for multimodal imaging studies on the same group of subjects, using the same methods 
of network construction, and employing the same graph diagnostics. Complementary 
theoretical work, using for example neural mass models [30l [61] [83] , could potentially 
explain mechanistically why one might expect different disease-related changes in network 
structure in different functional or structural imaging modalities. 

Dynamic Network Variability 

We explored the evolution of functional connectivity during the performance of a 
working memory task. This exploration was guided by theoretical work proposing the 
critical role of functional network reconfigurations in healthy cognitive function |31j and 
by empirical work demonstrating the importance of reconfiguration flexibility for cogni- 
tive phenomena like learning |16j . We identified greater temporal network variability on 
average in people with schizophrenia than in healthy controls. This finding is consistent 
with theoretical formulations of dysconnectivity in schizophrenia that suggest a temporal 
disorganization of sequentially expressed dynamical states [18] . Indeed, from a theoreti- 
cal point of view, the inability of the dynamic brain network to sustain the organization 
necessary for healthy cognitive function might empirically result in the fragmentation of 
cortical and subcortical networks seen in schizophrenia [87 . Empirically, our findings are 
also consistent with a recent EEG study by Schoen and colleagues [SU] demonstrating 
an increased entropy of connections in schizophrenia specifically in the high frequency 7 
band. Together, these results point toward altered temporal trajectories of whole-brain 
cortical function potentially due to inadequate or altered constraints on brain dynamics 
in schizophrenia. 

Methodological Limitations and Considerations 

In addition to intra-frequency investigations, we have constructed inter-frequency 
networks to determine whether cross-frequency communication patterns are altered in 
schizophrenia. Cross-frequency interactions are thought to facilitate the cross-function 
integration of information over spectrally distinct processing streams |65j and might be 
critical for memory function [79] . Phase-amplitude coupling and similar methods [251 151 [5] 
have been used to quantify these interactions, but the theoretical foundations of these 
methods are still not well understood [90] E2] ■ We have chosen to use a relatively simple 
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nonlinear measure of interaction - the mutual information between times series extracted 
from different wavelet bands - to quantify statistical associations between frequencies 
that could be studied from a network perspective. Although outside of the scope of 
the present paper, it would be interesting in future to examine the effects of alternative 
estimates of these interactions on network structure. 

The data used for this study was acquired using a CTF machine, whose axial gra- 
diometers have source profiles that include information from a wide spatial range, limiting 
the potential for anatomical interpretations. Commonly used source localization tech- 
niques allow for greater confidence in anatomical localization but simultaneously change 
the correlation structure between time series. Instead, to retain the inherent network 
correlation structure and increase the localization specificity, we transformed the data 
into planar space [T31 |S5] . 

In the present study, we did not observe significant correlations between the single- 
valued diagnostics (entropy, strength, cost-efficiency and variability) and the accuracy of 
task performance, for which our small sample size might be a factor. A more extensive 
examination of the relationship between behavioral variables and network diagnostics is 
outside of the scope of this study, whose focus was primarily to identify alterations in 
network structure and dynamics in schizophrenia. In this study, we have reported results 
at the statistical threshold of p < 0.05 uncorrected for multiple comparisons. Our choice 
is justified given the exploratory nature of the work, and is moreover supported by the 
fact that it is not yet known how to correct for comparisons across multiple network 
diagnostics. The interdependence of these measures is not well understood, but prelim- 
inary evidence suggests the presence of non-trivial correlations between them |57l 114] , 
such that strict multiple comparisons correction would not be appropriate. Furthermore, 
because all of the people with schizophrenia included in this study were on medication, 
most were smokers, and their performance on the task was significantly worse than the 
performance of the normal subjects, we can not determine whether our findings are 
driven by the disease, by associated epiphenomena, as an effect rather than a cause of 
the poor task performance (e.g. altered attention, effort, distraction) or by medication or 
a combination of these. To examine the effects of these various epiphenomena, it would 
be important to perform a follow-up study in unaffected siblings, for which our study 
gives important guiding information. 

Conclusion 

Recent advances in physics and mathematics have provided unique, robust quantita- 
tive network methods to examine the structure and organization of whole-brain functional 
connectivity. Mounting evidence from a plethora of imaging modalities, cognitive states, 
and diseases underscores the utility of network theory in capturing previously uninves- 
tigated variations in large-scale brain function and its alteration in disease states [§]■ 
However, the interpretation of these findings is complicated by the fact that network 
signatures of schizophrenia are not consistent across imaging modalities and analytic 
methods and only provide indirect information about the actual mechanisms responsible 
for the results. In this study, we have introduced several new methods of studying dy- 
namic properties of brain networks in schizophrenia, including cross-frequency networks 
and temporal variation in network structure, and we hope that they will contribute to 
such advancements. 
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Supplementary Material for Intra- and Inter-Frequency Brain Network Struc- 
ture in Health and Schizophrenia 

This supplementary materials section includes the following: 

• Mathematical definitions for network diagnostics employed in the present analysis 

• Supplementary Figure 1: Entropy and Strength 

• Supplementary Figure 2: Network Diagnostics: Part I 

• Supplementary Figure 3: Network Diagnostics: Part II 

• Supplementary Figure 4: Variability of Network Diagnostics: Part II 
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Mathematical Definitions 
Time Series Diagnostics 

Wavelet Entropy: Entropy, as denned by Shannon |82| is a simple, time-independent 
measure of the degree of order/disorder of the signal and has been applied to measure the 
univariate complexity of signals obtained in neuroimaging [731 114) . Here, we calculated 
wavelet entropy with the MATLAB function wentropy.m: 

£00 = -E s ? lo s( s ?) ( 2 ) 

i 

where s is the signal of a single region in a given individual and Sj are the coefficients of 
s in the orthonormal wavelet basis. 



Network Diagnostics 

A network is composed of units (nodes) and connections between those units (edges). 
The degree ki of node i is defined as the number of edges emanating from node i. 

Clustering coefficient The clustering coefficient C is defined by supposing that a node 
i has ki neighbors, so a maximum of ki(ki — l)/2 edges can exist between these neighbors 
[55] . The local clustering coefficient Cj is the fraction of these possible edges that actually 
exist: 

ki(ki — 1) 

The clustering coefficient C of an entire network is then defined as the mean of d over 
all nodes i. 

Hierarchy A sense of hierarchical structure of the network can be characterized by 
the coefficient /3, which is a parameter quantifying the putative power law relationship 
between the clustering coefficient Cj and the degree ki of all nodes in the network [71] : 

Ci-kf. (4) 

Pragmatically, we estimate j3 using the best linear fit of C vs. k in loglog space with a 
robust outlier correction. 

Assortativity The degree assortativity of a network (which is often called simply 
'assortativity') is defined as 

where ji and ki are the degrees of the nodes at the two ends of the i th edge (i 6 {1 , . . . , E} 
[53]. The assortativity measures the preference of a node to connect to other nodes of 
similar degree (leading to an assortative network, r > 0) or to other nodes of very different 
degree (leading to a disassortative network, r < 0). Social networks are commonly 
found to be assortative while networks such as the internet, World-Wide Web, protein 
interaction networks, food webs, and the neural network of C. elegans are disassortative. 

Mean Connection Distance The estimated connection distance of an edge, dij, is 
defined as the Euclidean distance between the centroids of the connected regions i and j 
in standard stereotactic space. The mean connection distance, d, is defined as the average 
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connection distance over all edges in a network [TU]. Thus connection distance differs 
from the other, topological and dimensionless network diagnostics in that it represents a 
spatial or topographic property of the network and has units of distance (mm) . 

Rent's exponent Rent's exponent is a topophysical property of a network; that is, it 
describes how a non-physical topology is embedded into a physical space, which in the 
case of neuronal fiber tracts is the physical space of the brain [TTJ. Rent's rule, which 
was first discovered in relation to computer chip design, defines a scaling relationship 
between the number of external signal connections (edges) e to a block of logic and the 
number of connected nodes n in the block [37J : 

e ~ n p , (6) 

where p G [0,1] is the Rent exponent. Following [TT], the Rent's exponent is found by 
tiling the Euclidean space of the network with Nb ox = 5000 overlapping randomly sized 
boxes (e.g., three-dimension cubes). In each box we determine the number of nodes 
(n) and the number of connections (e) that cross the box boundaries. The gradient 
of a straight line fitted to log(n) versus log(e) using iteratively weighted least squares 
regression is an estimate of the Rent exponent p. To minimize boundary effects, p is 
estimated using the subset of boxes which contains less than half the total number of 
nodes, n < N/2. 

Global efficiency The global efficiency was defined by Latora and Marchiori 53 and 
first applied to ncuroimaging data in [1 . The regional efficiency of a single node, i, is 
defined as 

^) = ^E^' (?) 

where i = 1, 2, 3, . . . , N indicates the index region, j ^ i denotes a region connected to 
i, and L^j is the minimum path length between regions i and j. Regional efficiency is 
therefore inversely related to minimum path length and a region with high efficiency will 
have short minimum path length to all other regions in the graph. The global efficiency 
of a graph is defined as the mean of E{i) over all possible regions, and is commonly 
denoted E g i i,. 

Local efficiency Latora and Marchiori also defined a local efficiency, which measures 
the efficiency of the subgraph surrounding node i: 

where Gi is the subgraph of nodes and edges connected to node i and Lj^ is the minimum 
path length between nodes j and k in the subgraph |53j . 

Betweenness centrality Geodesic node betweenness or more simply 'betweenness cen- 
trality' is defined for the i th node in a network Q as 



Bi= £ ?f^, (9) 

where all three nodes (j, m, and i) must be different from each other, V'j.m is the number 
of geodesic paths between nodes j and m, and ipj, m (i) is the number of geodesic paths 
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between j and m that pass through node i. The betweenness centrality of an entire 
network B is defined as the mean of Bi over all nodes i in the network. 

Modularity Networks can be partitioned into communities or modules 68, 3S] where 
nodes inside the same community are more densely connected to each other than they 
are to nodes in other communities. The modularity 43, 63 of a network partition is 
defined as: 

where fc, is the degree of node i, m is the total number of edges in the network, is an 
element of the adjacency matrix, 5ij is the Kronecker delta symbol, and Cj is the label 
of the community to which node i has been assigned |54j . Here we used the Louvain 
locally greedy algorithm [17] to optimize the modularity quality function over the space 
of possible network partitions. We report the maximum value of Q over this optimization 
procedure. 

Synchronizability The synchronizability, S, of a network characterizes structural prop- 
erties of a graph that hypothetically enable it to synchronize rapidly [T2]. The synchro- 
nizability is defined as 



where A2 is the second smallest eigenvalue of the Laplacian C of the adjacency matrix, 
and A jv is the largest eigenvalue of C [TJ. 

Robustness The robustness metric, p, indicates the network's resilience to either tar- 
geted, p t , or random, p r , attack. In a targeted attack, hubs are removed one by one in 
order of degree, k, while in a random attack, nodes are removed at random independent 
of their degree. Each time a node was removed from the network, we re-calculated the 
size of the largest connected component, s. Robustness is then usually visualized by a 
plot of the size of the largest connected component, s, versus the number of nodes re- 
moved, n [HH7J- The robustness parameter, p, is defined as the area under this s versus 
n curve. More robust networks retain a larger connected component even when several 
nodes have been knocked out, as represented by a larger area under the curve or higher 
values of p. 

Cost Efficiency We define the cost efficiency [TJ [T3J EZ] at a node as the maximal 
difference between the regional efficiency and the network density or cost over the inves- 
tigated range of cost values: 

CE(i) = E(i) - K. (12) 

We also define the cost efficiency of the network CE net as the mean of CE(i) over all 
nodes in the network. 
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Figure 9: Entropy and Strength. Correlation between the strength of connectivity and complexity, as 
measured by wavelet entropy, in 7-, ft-, a- and 6>-bands. Single data points represent the mean value pairs 
from a block of trials, trials from different subjects are distinguished by different colors and markers. 
Red, orange and pink markers denote subjects with schizophrenia spectrum diagnosis; blue, turquoise 
and purple markers denote healthy subjects. Value pairs from single subjects exhibit a tendency to 
appear in clusters, which (mostly) are broken up only at low entropy and/or high connectivity. With 
some subjects, value pairs appear outside the main cluster in all bands, with others, only in some. Also 
displayed are fitted linear functions for the two groups (red and blue lines) with r 2 values as indicators 
of goodness of fit. It should be noted that these were obtained by fitting to the mean value pairs for 
subjects, averaged over all 6 blocks of trials. This was done to avoid fitting to values for which there are 
two sources of variance (subjects and blocks). These fits indicate a negative correlation between entropy 
and strength, especially for healthy subjects where r 2 values are much higher. 
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Figure 10: Network Diagnostics: Part I One set of six network diagnostics (global efficiency, be- 
tweenness centrality, clustering coefficient, local efficiency, modularity, and hierarchy) is plotted as a 
function of density in networks within and between frequency bands. Each curve represents one subject, 
values averaged over all 66 trials. Curves for healthy controls are black, those for SZ patients colored. 
The two sets of curves were tested for statistically significant difference with Functional Data Analysis 
(FDA), the resulting p-values are given. Where significance (p < 0.05) was calculated, the color of the 
SZ curves was set to red, purple otherwise. We see significant differences for most diagnostics between 
the groups in the 7, /3 and a bands, as well as for the 7 — cross-frequency network. 
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Figure 11: Network Diagnostics: Part II A second set of six network diagnostics (synchronizability, 
assortativity, robustness to targeted and random attack, Rent's exponent, and mean connection distance) 
is plotted as a function of density in networks within and between frequency bands. Each curve represents 
one subject, values averaged over all 66 trials. Curves for healthy controls are black, those for SZ patients 
colored. The two sets of curves were tested for statistically significant difference with Functional Data 
Analysis (FDA), the resulting p-values are given. Where significance (p < 0.05) was calculated, the 
color of the SZ curves was set to red, purple otherwise. In the 7 and ft bands, as well as in the 7 — 
cross-frequency network, we see again a majority of diagnostics showing significant (p < 0.05) differences, 
but not in the a band. 
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Figure 12: Variability of network diagnostics. Coefficient of variation for binary network diagnostics 
in all intra- and inter-frequency networks. Values indicate variability over trials, averaged over all healthy 
(blue) and schizophrenic subjects (blue) and over the entire range of cost values. Error bars indicate the 
square mean of the standard errors over subjects and costs. 
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